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Abstract. It is well-known that the convergence of Krylov subspace methods to solve linear 
system depends on the spectrum of the coefficient matrix, moreover, it is widely accepted that 
for both symmetric and unsymmetric systems Krylov subspace methods will converge fast if the 
spectrum of the coefficient matrix is clustered. In this paper we investigate the spectrum of the 
system preconditioned by the deflation, coarse correction and adapted deflation preconditioners. 
Our analysis shows that the spectrum of the preconditioned system is highly impacted by the angle 
between the coarse space for the construction of the three preconditioners and the subspace spanned 
by the eigenvectors associated with the small eigenvalues of the coefficient matrix. Furthermore, 
we prove that the accuracy of the inverse of projection matrix also impacts the spectrum of the 
» . preconditioned system. Numerical experiments emphasized the theoretical analysis. 
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1. Introduction. We consider the iterative solution of a linear system 

Ax = b, 



where A 6 R" x ™ is symmetric and positive definite (SPD). It is well-known that 
the convergence of Krylov subspace methods for solving linear systems depend on 
the eigenvalue distribution of A. Recently, several studies [1 [3 H ITUrfiil [13] 
have shown that by removing the subspace spanned by the eigenvectors corresponding 
to several small eigenvalues from the Krylov search space makes the spectrum more 
,__! clustered and consequently the convergence is improved. 

^ In this paper, we refer to matrix of the form 

o 

^H (1.1) E = Z T AZ, z e R nxr . 

o 

as a projection matrix, and the subspace spanned by the columns of Z is referred 
^f to as a coarse space. In an ideal situation, the coarse space contains the vectors 

corresponding to the lower part of the spectrum that is responsible for the stagnation 

of Krylov subspace methods. As shown in this paper, the preconditioned systems will 

have the desired properties when the preconditioner is enchanced with the ideal coarse 

^ subspace. In contrast to the general coarse space, we refer to the coarse space spanned 

l> by the eigenvectors associated with several small eigenvalues of A as the exact coarse 

space. 

Next we briefly mention a few existing approaches that take the form of a standard 

precondition with an additive coarse space enchancement. In [21 E2], the deflation 

preconditioner is defined by 

(1.2) P D =1- AZE- X Z T . 

Obviously PdA is singular since Pd is singular. Fortunately, Krylov subspace methods 
converge for singular linear systems as long as they are consistent; furthermore, zero 
eigenvalues do not impact the convergence since the corresponding eigenvectors never 
enter the Krylov subspace [T2] . 
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Instead of zero out the small eigenspace in the deflation method, the precondi- 
tioned based on coarse correction shift the small eigenvalues to the large ones. In 
[19) , the coarse correction preconditioner in domain decomposition method is defined 

by 

(1.3) P c = I + ZE- l Z T . 

The abstract additive coarse correction is M~ l + ZE~ l Z T , where M is the sum 
of the local solves in each subdomain. Another popular preconditioner in domain 
decomposition method is the abstract balancing preconditioner [3] 

(1.4) Pbnn = (I - ZE^ 1 Z T A)Ar 1 {I - AZE^ 1 Z T ) + ZE^ 1 Z T . 
The adapted deflation preconditioner [TS] is defined as 

Padefi = M- l P D + ZE- X Z T . 

It is shown in |19) that Padefi is cheaper than Pbnn but is as robust as Pbnn- 
Moreover, PbnnA and PadefiA have an identical spectrum. In both Pbnn and 
Padefi, the first term corresponds to the fine space and the second term is for the 
coarse space, therefore they are called two-level preconditioners. Throughout this 
paper, we restrict our analysis to one-level methods. Let M in Padefi be /. We 
define Pa as 

(1.5) P A =1- AZE^Z? + ZE- l Z T . 

Let X be any basis of a coarse space not Z. Obviously the following identity 

X(X T AX)- 1 X T = z{z T Azy l z T 

holds in exact arithmetic since there exists a nonsingular matrix C £ M. rxr such that 
X = ZC. This implies that Pd, Pc and Pa are determined uniquely by the coarse 
space. Thus we can choose an appropriate basis to form Z and then to construct a 
preconditioner with certain desirable properties. 

As we will see, if an approximate coarse space is used to construct precondition- 
ers, then the spectrum of the preconditioned systems is related to the angle between 
the approximate and exact coarse spaces. We also prove that the coarse correction 
and adapted deflation preconditioners are more roust than the deflation precondi- 
tioner when the projection matrix is solved inexactly. In section 2, we first review 
the spectral properties of the preconditioned system when the preconditioners are 
constructed with the exact coarse space, then we estimate the spectral bounds of the 
preconditioned systems when the approximate coarse space is used. Section 3 presents 
the perturbation analysis on the spectrum of the preconditioned system in the case 
that the projection matrix have some perturbation. Numerical results are reported 
in Section 4. 

2. Coarse space spanned by the approximate coarse space. In this sec- 
tion, we briefly review the spectrum of the system preconditioned by using the exact 
coarse space, then based on these properties we derive the bounds of the spectrum of 
the system preconditioned by using the approximate coarse space. Let (Aj,«j) be an 
eigenpair of A and vt be normalized, (vi, ■ ■ ■ ,V n ) is orthogonal since A is SPD. The 
spectral decomposition of A can be written as 

<"> *-<™>(o£)(? 



where A = diag{Ai, • • • , A r }, Aj_ = diag{\ r +i, • • • , A n }, V — (v\, ■ ■ ■ ,v r ), and V± = 
Assume that Ai, A2, ■ ••, A,- are small eigenvalues that impacts the convergence 



of the Krylov subspace methods. Theorem 2.1 and Theorem 2.2 show that if we take 



Z = V, the small eigenvalues are removed or shifted when the preconditioners are 
applied on either side of A, moreover, the rest of eigenvalues and all the eigenvectors 
are not changed. 

Theorem 2.1. Let E = V T AV. Define 

P D =I- AVE^V 7 \ 

P C = I + VE- 1 V T , 

P A = I - AVE^V 7 + VE~ 1 V T . 

Then we have the following spectral decomposition 

^-<™>(U)(S 

Proof. From the definition of E, we have 

E = V T AV = V T VA = A. 

Next, let us consider the first result. Obviously, PdAV = 0. We then have 

P D AV± = AV ± - AVA- 1 V T AV ± 
= Vj_A± - VV T Vj_A ± 
= Vj_A ± . 

Thus the first spectral decomposition holds. Since 

P C A = A + AVE- l V T = A + ^AA" 1 ^ = A + VV T , 

we have 

P C AV = AV + VV T V = V(A + I) 

and 

P C AV± = AV X + W T V X = AV X = V x Xx- 

Hence, the second spectral decomposition is true. The third spectral decomposition 
follows from 

P A AV =AV- AVE- 1 V T AV + VE^V T AV = V 

and 

P A AV± = AV ± - AVE- 1 V T AV ± + VE- 1 V T AV ± 

= AV± - AVE~ 1 V T V A _A A _ + VE- 1 V T V A _A A _ 
= AV± = V± A j. 



It follows immediately from these three spectral decomposition that PdA, PcA and 
Pa A are symmetric. D 



Theorem 2.2. Pp, Pc and Pa are defined in Theorem 2.1 Then we have 
AP D = PdA, AP c = P C A and AP A = PaA. 
Proof. The first result follows from 

AP D V = AV- AVAA~ 1 V T V = = P D AV 

and 

AP d V_l = AV X + AVA~ 1 V T V 1 _ = V X A± = P D AV ± . 

Note that A, Pc and PcA are symmetric. We then have 

P C A = (P c Af = A T pr = AP C . 

Since 

AP A = A - A 2 VE- 1 V T + AVE^V 7 

= A - AVAA - 1 V T + VAA~ 1 V T 
= A - AVV T + VV T , 

we get 

AP A V = AV - AVV T V + VV T V = V= P A AV 

and 

AP A V± = AV± + AVV T V±_ + VV T V±_ = AV±_ = V±A ± = P A AVj_. 

The third result is therefore true. In addition, APjj , APc and APa are symmetric as 
well. D 

For a large system, it is impractical to build the preconditioners by using the 
exact eigenvectors associated with the small eigenvalues, since in general computing 
these eigenvectors is more costly and more difficult than solving a linear system. 
Nevertheless, for some cases, the approximate eigenvectors are cheaply obtained and 
thus are used to produce the preconditioners. For example, in the Newton method for 
solving nonlinear problems, the information obtained during solving the first linear 
system is reused to build the coarse spaces for accelerating the convergence of the 
succeeding linear systems, see [7] and [5] . We hope that the system preconditioned by 
using the approximate coarse space has the similar eigenvalue distribution as described 
in Theorem |2.1| This motivates us to analyse how the perturbation in the coarse space 
impacts the spectrum of the system preconditioned by Pd, Pc and Pa- 

In the following discussion we assume that Z is column orthogonal. Let Z be the 
subspace spanned by the columns of Z . Let Z 1 - be the orthogonal complement of Z 
and Z± be an orthogonal basis of Z 1 - . Likewise, let V be the subspace spanned by the 
columns of V, V ± be the orthogonal complement of V and V± be an orthogonal basis 
of V ± . Let a denote the singular value of a matrix. Let dist(Z,V) denote the distance 
between subspaces Z and V. It is shown in jS] that dist(Z,V) can be evaluated by 
either a max (Z T V±) or a max (V T Z±). Note that < dist(Z,V) < 1 since Z and V 
have the same dimension. We define the acute angle between subspaces Z and V as 



9 = arcsin dist(Z, V). The next lemma shows that cos9 can be evaluated in a similar 
way as sin 9. 

Lemma 2.3. Let 9 be the acute angle between subspaces Z and V that have the 
same dimension. Let Z and V be the orthogonal bases of Z and V respectively. Then 

sine* = o- max (Z T V±) = o- max (V 7 , Z±), 
cos6 = o- mm (Z T V) = o- mm (Z 1 [V±). 

Proof. The first identity and its proof can be found in [6l Theorem 2.6.1]. We 
only prove the second identity here. Since (V, VI) is orthogonal, it follows from 
|| (V, V±) T Zx\\ 2 = 1 for all unit 2-norm x G R r that \\V T Zx\\l + \\V[Zx\\l = 1. Thus 

o- m m(V T Z) 2 = min || V T Zx\\ I = 1 - max \\V[Zx\\j 

11*112=1 ||*||2=1 

= 1 - O-maAVl Zf = COS 2 9. 

Similarly, since (Z, Z±) is orthogonal, it follows from \\(Z, Z±) T V±x\\ 2 — 1 for all unit 
2-norm x G R^"^ that \\Z T V_ L x\\ 2 + \\ZlV x x\\% = 1. Thus 

<? min (ZlV ± ) 2 = min ||ZjV x a;||| = 1 - max \\Z T V ± x\\ 2 

11*112 = 1 ||a;||2 = l 

= l-a max (Z T V ± ) 2 = cos 2 9. 
The second identity is therefore valid. D 



Let Pjj be defined by ( 1.2 1. PjjA is an symmetric matrix since A is SPD. Hence, 
Courant-Fischer Minimax Theorem [5] can be applied to estimate the eigenvalues of 
PdA. As is known that if A G R™ xn is symmetric, then 



Xk{A) = max min x Ax, k = 1, 

dim(S) — k x£»S,||a:||2=l 



Theorem 2.4. Let Pd be defined by (1.2) and 9 be the acute angle between V 
and Z. Then PdA has r zero eigenvalues. Moreover, if cos 9 ^ 0, then the nonzero 
eigenvalues of PdA satisfy 

A mi „(A_L) - e D < X(PdA) < A max (Ax) + vd, 

whererjD = A maa; (A_ L )(sin6'+sin 2 9) ande D = r] D + \\E^ 1 \\ 2 (\\E\\ 2 +\ m ax(A±)) 2 tan 2 9. 
Proof. It follows from PdAZ = that PdA has r zero eigenvalues. From Theorem 
|2.1| we have 

(2.2) x T P D Ax = x T V±A±V[ x. 

For all unit 2-norm x G R n , it can be written as x = X\ + x 2 , where xi G Z 
and x 2 G Z- 1 . Moreover, there exist t G M. r and s G M. n ~ r such that x\ = Zt and 
x 2 = Z±s, and 

x T AZE- 1 Z T Ax = XiAZE~ 1 Z T Axi + x^AZE- 1 Z T Ax 2 
+x%AZE- 1 Z T Ax 1 +x 2 r AZE- 1 Z T Ax 2 
= t T Z T AZt + t T Z T Ax 2 + x^AZt 

+x^AZE- 1 Z T Ax 2 
= x\ Ax\ + x[Ax 2 + x 2 Axi + x 2 AZE~ 1 Z T Ax 2 
= x T Ax - x 2 Ax 2 + x 2 AZE~ 1 Z T Ax 2 . 



Then we obtain 
(2-3) 



x T P D Ax = x T Ax - x T AZE~ 1 Z T Ax 

= xjAx 2 - x^AZE~ 1 Z T Ax 2 . 



Subtract (2.2) from (2.3) on both sides, we get 

(2.4) x T P D Ax = x T P D Ax + x 2 r Ax 2 - x T V±AxVfx - x\ AZE^ 1 Z T Ax 2 . 

The middle terms on the right-hand side of the above expression can be replaced by 

x 2 Ax 2 - x T V±_A±y1x = x 2 Ax 2 - (xj + x 2 )V±A±Vl(x 1 + x 2 ) 
= x^VAV T x 2 - x\ VlAxVf xi 

-xJV ± A ± V[x 2 - x%V±A 1 ylxi 
= s T (ZlV)A(V T Z ± ) s - t T (Z T V ± )A ± (V?Z)t 

-t T (Z T V 1 _)A^{VlZ 1 :)s - s T {ZlV x )A 1 _{VlZ)t. 

Il^ilb = Plb and HX2II2 = ||s||2 since Z and Z±_ are column orthogonal. Using Lemma 
|2.3[ we obtain 

\xlAx 2 -x T V x A x Vlx\<{\\x 2 f 2 \\A\\ 2 + \\x^ 2 \\A x \\ 2 )^B 

+2 1 i ^! 1 1 2 1 1 ^ 2 1 1 2 1 1 A^ 1 1 2 sin 6» cos 6>. 



(2.5) 



For the last item on the right-hand side of (2.4), we only need to estimate the 
bound of Z\AZ because of x^AZE' 1 Z T Ax 2 = s r (Z\AZ)E- x {Z\AZ) T s. It follows 
from (Z, Z X ){Z, Z ± ) T = I that 

AZ = Z(Z T AZ) + Z±(Zj_AZ). 

Multiply by Vf from the left on both sides of the above equation, we have 

(V[Z ± )(ZlAZ) = V[AZ - {V[Z)(Z T AZ) 

= A ± {V[Z) - {V[Z)(Z T AZ). 

Vj_Z\_ is invertible since cos 9 7^ 0. Thus 

(2.6) (ZlAZ) = {VlZ^A^VlZ) - (V[Z 1 _)- 1 (V[Z)E. 



Using Lemma 2.3 we have \\Z]_AZ\\ 2 < (\\E\\ 2 + ||A_|_|| 2 ) tan0. Hence 

(2.7) x^AZE- 1 Z T Ax 2 < \\x 2 \\l\\E-%(\\E\\ 2 + ||A_l|| 2 ) 2 tan 2 0. 
E^ 1 is SPD since A is SPD. Consequently, 

(2.8) x 2 AZE~ 1 Z T Ax 2 = (Z T Ax 2 ) T E- 1 (Z T Ax 2 ) > 0. 

PdA is symmetric since A is SPD. Applying Courant-Fischer Minimax Theorem 
to ( |2~4"| ) with ( |2~5| ) and ( [278] ), we have 

\{P D A) < X(P D A) + \x T 2 Ax 2 - x T V±A 1 Vlx\ 

< A moa; (Ax) + A moa: (Aj.)(sin0 + sin 2 0). 



Note that PjjA has r zero eigenvalues. Again, applying Courant-Fischer Minimax 



Theorem to (2.4) with (2.5 1 and (2.7), the lower bound of the nonzero eigenvalues of 



PdA is given by 

X{P D A) > X(P D A) - \xlAx 2 - x T VxkiVl x\ - x%AZE- 1 Z T Ax 2 
> A mm (A ± ) - \ max (A ± ) (sin 6 + sin 2 9) 

-\\E- 1 \\ 2 (\\E\\ 2 + x max (A ± )) 2 tan 2 9. 

As a result the theorem is true. D 

The above theorem shows that in exact arithmetic, as 9 approaches zero, the max- 
imal and minimal nonzero eigenvalues of PdA converge to X max (A±) and A m i„(Aj_) , 
respectively. Hence with an appropriate coarse space the spectrum of PdA is similar 
to that of PdA. When there exists rounding error, however, PdAZ may not be equal 
to a zero matrix. In this case, PdA possibly has some eigenvalues around zero that 
should be equal to zero in exact arithmetic. So there is a potential risk for Pd to 
yield a poor spectrum of the preconditioned system. 

The authors in [T!5] investigated the properties of Pd, Pbnn and Padefi- They 
established the relations between these preconditioners in terms of the spectrum. 
Suppose that M is an SPD matrix. Let the spectrum of PdM^ x A be given by 
{0, . . . , 0, 7 r+ i, . . . , 7„} with 7 r+ i < 7,. +2 < • • • < 7„. Let the spectrum of P B nnA 
and PadefiA be {1, . .. , 1,/V+i, • • • ,Mn} with ^ r+1 < /i r+2 < • • • < /J, n . Then, 
7i = Hi for all i = r + 1, . . . ,n. The proof of this result can be found in [THl Theorem 
3.3]. From the relation of Pd and Padefi-, we immediately obtain the next corollary 
if we take M = I. 

Corollary 2.5. Let the spectrum of PdA be given by {0, . . . , 0, A r+ i, . . . , A„}. 
Then the spectrum of PaA is {1, . . . , 1, A r+1 , . . . , A„}. 

Corollary |2 . 5 1 implies that if the spectrum of PdA is known, the one of Pa A would 



be known. The spectral bounds of PdA are described in Theorem 2.4 so we can easily 
bound the spectrum of PaA. 



Theorem 2.6. Let Pa be defined by {1.5). Let 9 the acute angle between subspaces 



Z and V. If cos 9 ^ 0, then the eigenvalues of PaA satisfy 

min{l, A m j„(Aj_) - e D } < X(P A A) < max{l, X max (A ± ) + n D }, 



where r\D and £d ar ^ defined in Theorem \2.4\ 

Proof. It follows directly from Theorem |2.4| and Corollary |2.5| D 
Next, we consider the spectrum of Pc'A. Note that PqA is not necessarily sym- 
metric although both A and Pc are symmetric. So we can not apply Courant-Fischer 
Minimax theorem to estimate the eigenvalues of PqA, which are positive since PqA 
is similar to a SPD matrix, and are a subset of x T PqAx for all unit 2- norm x £ M. n . 
Hence the eigenvalues of PqA can be bounded by estimating x T PqAx. 



Theorem 2.7. Let Pc be defined by (1.3), and 9 the acute angle between V and 
Z. If cos 9 7^ 0, then 

Xmax(PcA) < max{l + A mox (A), A mox (Aj_)} + £ C] 
X m in{Pc'A) > min{l + X min (A) , X m in(A±)} - e c , 



2 , 



where ec — \{X ma x{A±)\\E 1 || 2 + 1) tan# + sin# + sin 

Proof. P C A is similar to A + A 1 / 2 ZE^ 1 Z T A 1 / 2 since A is SPD. Moreover A + 
A 1 / 2 ZE- 1 Z T A 1 / 2 is SPD as well since A and E are SPD. Thus the eigenvalues of 
PcA are positive. 



For all unit 2- norm x € K n , we write x — X\ + x%, where X\ € Z and x 2 G Z±. 
There exists f gK r such that xi = Zt, likewise, there is s £ M. n ~ r such that x 2 = Z±s. 
Then x T PcAx can be expressed as follows 



(2.9) 



x T P c Ax = x T Ax + (xi + x 2 ) T ZE~ 1 Z T A(x 1 + x 2 ) 
= x T Ax + x^ZE- 1 Z T Ax 1 + x\ZE- x Z T 
= x T Ax + x[x x + x[ ZE- 1 Z T Ax 2 - 



From the definition of Pc in \2.1\, we obtain 
(2.10) 



x T P c Ax = x T Ax + x T VE- 1 V T Ax 



= x T Ax + x T VV T x. 



Subtract ( |2.10| from ( |2.9| ), we have 

(2.11) x T P c Ax - x T P c Ax = xJZE- 1 Z T Ax 2 + x{xi - x T VV T x 

= xJZE- 1 Z T Ax 2 + x[VxV[xi 



-x 2 VV T x 2 -2x r [VV T x 2 . 



Since both A and E arc symmetric, 



xl ZE~ 1 Z T Ax 2 = x 2 AZE- 1 Z T x 1 = s^Z^AZE^t. 



Using (2.6), we have 

(ZlAZ)E-i = (VlZ 1 _r 1 Ax(VlZ)E- 1 (V[ Z i)" 1 (V[ Z) . 



Note that ||aJi||2 = ||i||2 and ||#2||2 = ||s||2- Using Lemma 2.3 we obtain 
(2.12) \xlZE- l Z T Ax 2 \ < ||a;i||2||a;2||2(||Aj.||2||^- 1 ||2 + l)tan^ 

Since ||x||2 = 1, we have the following bounds with Lemma 2.3 



(2.13) 



Q<XiV±V[xi < ||xi||^sin 2 6l, 

0<x%VV T x 2 < ||a;2||l sin 2 6», 

\x{VV T x 2 \< ||si|| 2 ||x2||2 sin0. 



With ( |2.11[ ), ( |2.12| and ( |2.13[ ), we obtain 

\ min {PcA) -e c < x T P c Ax < \max{PcA) + e c , 

where ec — j(A maa ;(Aj_)||£ ,_1 ||2 + 1) tan 8 + sin 8 + sin 2 8. Thus the theorem follows 
from mhx{x PcAx} < X(PcA) < max{a; T Pc Ax} for all unit 2-norm x € R™. D 

From Theorem |2.7| we conclude that the maximal and minimal eigenvalues of 
PcA converge to max{A ma:E (AjJ, 1 + A ma:E (A)} and min{A mm (Aj_), 1 + A min (A)} 
respectively as 8 approaches zero. With an appropriate coarse space, the spectral 
distribution of PcA would be close to that of PcA. Analogously, it can be proved 
that the maximal and minimal eigenvalues of APc have the same bounds as that of 
PcA. 



3. Inexact inverse of projection matrix. In practice, we do not form E^ 1 
explicitly to compute y = Ex. Here y and x are vectors with the suitable size. 
Instead, we compute the LU factorization of E once, then solve the two triangular 
linear systems to obtain y. But it is expensive to compute the LU factorization when 
the matrix E is large. We therefore replace E by a perturbed one that is cheaper to 
compute. In this seciton, we analyse how the perturbation in the projection matrix 
impacts the spectrum of the p reco nditioned matrix. Assume H is an invertible matrix 



as an approximation to E in (2.1 1. 

Theorem 3.1. Let P D = I - AVH~ 1 V T and p = EH^ 1 —I. If the eigenvalues 
of PdA are real, then 

-£d < HPdA) < X max (A ± ) + £ D , 

where £ D = \\p\\ 2 \\A\\2. 

Proof. \\p\\2 measures how much H is approximate to E when H~ l is used as the 
right inverse of E. Obviously ||p||2 = i f and only if E = H . Assume \\x\\2 = 1 and 



use the definition of Pr, in Theorem 2.1 we obtain 



x T P D Ax = x T Ax - x T AVH' l V T Ax 

= x T Ax - x T VV T Ax - x T VpV T Ax 
= x T P D Ax - x T VpAV T x. 

Since \x T V P AV T x\ < ||p|| 2 ||A||2, 

-HHhllAlb < x T P D Ax < \ max {P D A) + ||p|| 2 ||A|| 2 . 

Since the eigenvalues of PdA are real, the theorem follows from mm{x T Pd Ax} < 
X(PdA) < max\x T P D Ax}. □_ 



Theorem 3.1 states that PdA might have small eigenvalues around zero if ||p||2 ^ 
0, which leads to the worse spectral distribution than that of PdA. 

Theorem 3.2. Let P c = I + VH^V 1 ^ and p = H~ 1 E-I. If the eigenvalues 
of PcA are real, then 

Xnax{PcA) < max{l + A ma x(A), A max (AjJ} + £ c , 
^min(PcA) > min{l + X min (A), A mm (Aj_)} - £ c , 

where & = II p\\ i- 

Proof. In this theorem, we use H -1 as the left inverse of E. Assume ||x||2 
and use the definition of Pc in Theorem 



2.1 



then we have 



x T P c Ax = x T Ax + x T VH~ 1 V T Ax 

= x T Ax + x T VE~ 1 V T Ax + x T VpV T x 
= x T P c Ax + x T VpV T x. 

Since \x T VpV T x\ < \\p\\ 2 , 

Xmin(PcA) - \\p\\ 2 < X T P C Ax < Xmax(PcA) + \\p\\ 2 . 

Since the eigenvalues of PcA are real, the theorem follows from mm{x T PcAx} < 
KPcA) < maxjx T P c Ax}. U 



Theorem 3.2 implies that if the eigenvalues of PcA are real, the spectral distribu- 



tion of PcA is a little influenced with small ||p||2, because the maximal and minimal 



10 

eigenvalues of P C A converge to max{\ max (A_ L ),l + \ max (A)} and min{A mm (AjJ, 1 + 
A m in(A.)} respectively as p approaches zero. 

In next theorem, we need to estimate the difference of H and E when H^ 1 is 
used as both the left and the right inverse of E. 

Theorem 3.3. Let P A = I - AVh\- 1 V T + VH~ 1 V T . Let Pl = EH- 1 - I and 
pi = H~ l E — I . If the eigenvalues of Pa A are real, then 

\nax{PAA) < max{l,A moa .(A.]_)} + fA) 
X m in(PAA) > min{l, A m i„(Aj.)} - £4, 

«iAcret4 = ||pi||a||A||2 + ||p2||a- 

Proof. Assume ||x|| 2 = 1 and use the definition of Pa in (2.1), we have 

x T P A Ax = x T Ax - x T AVH- 1 V T Ax + x T VH- 1 V T Ax 
= x T Ax - x T V Pl AV T x + x T Vp 2 V T x 
-x T AVE~ 1 V T Ax + x T VE~ 1 V T Ax 
= x T P A Ax - x T V Pl AV T x + x T Vp 2 V T x. 

Since \x T V Pl AV T x - x T Vp 2 V T x\ < |M| 2 ||A|| 2 + ||p 2 || 2 , 

X mm (P A A) - ||pi|| 2 ||A|| 2 - ||p 2 || 2 < x T P A Ax < X max (P A A) + ||pi||a||A|| a + IHIa- 

We assumed that the eigenvalues of Pa A are real. Thus the theorem follows from 
mm{x T P A Ax} < X(PaA) < max{i T P A Ax} . D 

4. Numerical experiments. In this section, a numerical comparison of vari- 
ous preconditioners is reported. All tests are performed with Matlab (R2010b) on an 
Intel Core2 Duo E7500, 2.93GHz processor with 4Gb memory. Except for the defla- 
tion preconditioner, the system preconditioned by the coarse correction and adapted 
deflation preconditioners are not necessarily symmetric although A is SPD. Therefore 
we apply GMRES [TB] to solve the preconditioned system iteratively. In addition, we 
apply Gram-Schmidt method with reorthogonalization to maintain the orthogonality 
of basis of the Krylov subspace [6l [16] . 

4.1. Diagonal matrix. The first test case is a diagonal matrix with entries 10 -7 , 
1(T 6 , • • •, HT 1 , 1, 10, 10.1, 10.2, • • •, 209, 209.1. The matrix has 7 small eigenvalues 
less than 1 to be removed. The eigenvectors associated with these eigenvalues are 
the unit vectors, i.e., V — (e%, e 2 , . . . , e-?), where e^ is the zth column of the identity 
matrix. The right-hand side is a vector of all ones. The initial guess vector is a zero 
vector. All tests are required to reduce the relative residual below 10 -12 . GMRES 
method without preconditioning converges at 273th iteration. The perturbations in 
the coarse space and the projection matrix are generated by the Matlab function rand. 

Table |4~1~1 shows the distance between the exact coarse space and the coarse space 
with various perturbations. It should be noted that it is expensive and unnecessary 
in practice to compute sin 9 by Lemma |2.3| for a general linear system. Assume that 
(Aj, i>i) (i = 1, • • • , r) are Ritz pairs that are extracted from the perturbed coarse space 
by the Rayleigh-Ritz procedure [18 . In general, maxi = i i ... ir {||A«i — Aji>j||2} denoted 



by res max in Table 4.1 decreases as the two subspaces approach each other. So we 
can use res max to measure the distance between the two subspaces since it is more 
convenient to compute. 

In Table |4.2| the second column shows the number of GMRES iterations with 
the three preconditioners in the case that there is no perturbation in the coarse space 
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Table 4.1 
The distance between span{V} and span{V + rand/e}. 





e= le + 01 


e = le + 02 


e = le + 03 


e = le + 04 


£ = le + 05 


sin 8 


9.77e-01 


5.06e-01 


6.03e-02 


6.05e-03 


6.02e-04 


iCS-max 


7.08o+01 


5.54e+01 


7.33 


5.78e-01 


4.43e-02 



and the projection matrix. The columns 3-7 show the number of GMRES iterations 
when only the coarse space has some perturbation. As is shown, all preconditioners 
suffer from the perturbation if it is large (see the third column). As the perturbation 
decreases, Pc and Pa become better, whereas Pjj becomes better only when the 
perturbation is very small. On the other hand, if the perturbed coarse space is close 
enough to the exact one (see the last two columns), Pd is slightly more efficient than 
Pa, and both of them are more efficient than Pjj. 

Table 4.2 
The number of GMRES iterations with various preconditioners that are constructed with Z = 
V + rand/e and E~ 1 . 





V, E- 1 


£ = le + 01 


£ = le + 02 


£ = le + 03 


£ = le + 04 


£ = le + 05 


Pd 


71 


>300 


>300 


>300 


109 


88 


Pc 


104 


273 


267 


222 


173 


144 


Pa 


72 


290 


231 


167 


117 


96 



Figure [471] and Figure 4.2 show the eigenvalue distribution of the system precon- 
ditioned by the three preconditioners. As discussed in Section [51 because of rounding 
error, PjyA has some tiny eigenvalues around zero. In general, it is difficult to figure 
out the condition under which these tiny eigenvalues cause the stagnation in the con- 
vergence. For the test case of diagonal matrix, rounding error does not impact the 
convergence of PdA when the perturbation in the coarse space is significantly small. 
We also see that the spectrum of PaA is more clustered than that of PcA, which is 
consistent with the estimated bounds of their spectrum described in Theorem |2.6| and 
Theorem [271 
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Fig. 4.1. The eigenvalue distribution of the Fig. 4.2. The eigenvalue distribution of the 

preconditioned system. The preconditioners Pr), preconditioned system. The preconditioners Prj, 

Pc and Pa are built with Z = V + rand/le + 03 Pc and Pa are built with Z = V + rand/le + 05 

and E~ x . and E _1 . 



In Table [4~3| the first three rows show the number of GMRES iterations when only 
the projection matrix has perturbation. The difference of H^ 1 and E^ 1 is reported 
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in the last two rows. In comparison with the second column in Table |4~2j we conclude 
the perturbation in the projection matrix has a little impact on Pq and Pa, but has 
a severe impact on Po even when H~ l is very close to E _1 (see the last column). 

Table 4.3 
The number of GMRES iterations with different preconditioned, where only the projection 
matrix is perturbed and the perturbation of E is H = E + rand/e. 





e = le + 10 


e= le + 12 


e= le + 14 


e= le+16 


Pd 


>300 


>300 


>300 


>300 


Pg 


111 


104 


104 


104 


Pa 


88 


88 


80 


80 


\\H- X E - 1\\ 2 


1.68e-03 


1.08e-05 


1.77e-07 


1.23e-09 


WEH- 1 -!^ 


1.68e-03 


1.08e-05 


1.77e-07 


1.23e-09 



Figure [473] and Figure [474] report the spectral distribution of Pd, Pc and Pa- We 
see that Pc and P4 successfully shift the small eigenvalues of A to around one. Due 
to the perturbation in the projection matrix, Pd fails to deflate the small eigenvalues 
and thus yields some tiny eigenvalues around zero, which leads to a worse convergence 



(see the first row in Table 4.3 1 
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Fig. 4.3. The eigenvalue distribution of the Fig. 4.4. The eigenvalue distribution of the 

preconditioned system, where the preconditioners preconditioned system, where the preconditioners 

Prj , Pc and Pa are built with V and H = E + Pd , Pc ana ' Pa a-re built with V and H = E + 

rand/ le + 12. rand/le + 16. 



4.2. Boundary value problem. We solve the following model problem 

-V-(kVu) =/ in n= [0,1] 2 , 
u = on d£l, 

by the two-level multiplicative Schwarz method [8]. Here, k is the diffusion function 
of x and y. The model problem is discretized by FreeFem++ [T5] and the resulting 
coefficient matrix is of size 10201. Tests are performed on irregular overlapping de- 
compositions with the overlap of 2 elements. These overlapping decompositions are 
built by adding the immediate neighboring vertices to non-overlapping subdomain 
obtained by Metis [9]. 

In the two-level multiplicative Schwarz method, the first level preconditioner is 
the restricted additive Schwarz preconditioner (RAS) [T] that is responsible to remove 
high frequency modes of the original system, and the deflation, coarse correction and 
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adapted deflation preconditioners are applied as the second level preconditioners that 
remove lower frequency ones of the system preconditioned by one-level preconditioner 

inn. 

We choose Ritz vectors to span the coarse space, which are extracted from the 
Krylov subspace during the solve of the system preconditioned by RAS. These vectors 
are the approximate eigenvectors corresponding to the lower part of the spectrum of 
the preconditioned system. To enrich the information on lower part of the spectrum, 
we construct the coarse space by splitting Ritz-vectors, see [H [TTJ [19] and references 
therein. More precisely, let 



V = 



i'll 


«12 


■ Vl >T 


''21 


«22 


■ V 2 ,r 


unparts,! 


^nparts,2 


^npc 



store Ritz vectors columnwise, where nparts is the number of subdomains and r the 
number of Ritz vectors; let Zj store the orthogonal vectors obtained by orthogonalizing 
(Vji, Vi2, ■ ■ ■ , Vi jr ). Then Z is formed as follows 



Z = 



Z x 
z 2 







J nparts 



Obviously, span{T^} is a subspace of span{Z}; span{Z} is nparts times as large 
as spanjV^}. So span{Z} might have richer information corresponding to small eigen- 
values. 

Some comparisons of the three preconditioners are performed on two different 
configurations with highly heterogeneous coefficient k. See [TT] for details. Two cases 
are described as following: 

• skyscraper n: for x and y such that for [9x]=0(mod_2) and [9y]=0(mod 2), 
K = 10 4 ([9y] + 1); and k = 1 elsewhere. See Figure 
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• continuous k: n(x,y) = 10 6 /3sin(47r(x + y) + 0.1). Sec Figure 
For both cases, we use the zero vector as the initial guess vector, 
will stop when the relative residual is less than 10~ 10 . Moreover 
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The iteration 
we construct the 

coarse space with all the eigenvectors associated with the eigenvalues less than 0.5 
against the various domain decompositions. 

In Figures |4.7||4.10[ we report the convergence curves for the skycraper case 
against the various number of subdomains. Compared to one-level method, all three 
two-level methods improve convergence sufficiently. Pp and Pa have almost the same 
number of iterations although the initial residual of Pjj is much less than that of 
Pa- Pd and Pa are more efficient than Pq- Two-level method varies slightly on the 
number of iterations as the number of subdomains increases, while one-level method 
does. 

Figures 4.11|4.14 plot the convergence curves for the continuous case against the 
various number of subdomains. Once again, we see that two-level method with RAS 
and the three preconditioners all outperform one- level method with only RAS. Like 
the skycraper case, Pd and Pa have almost the same number of iterations for all four 
decompositions and outperform Pq- 
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Fig. 4.5. Skyscraper case 



Fig. 4.6. Continuous case 




Iteration number 



80 100 120 

Iteration number 



Fig. 4.7. Skyscraper case with 16 subdo- Fig. 4.8. Skyscraper case with 32 subdo- 

mains. 16 Ritz vectors spanning the coarse space, mains. 16 Ritz vectors spanning the coarse space. 
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Iteration number 



20 40 60 80 100 120 140 160 180 200 220 
Iteration number 



Fig. 4.9. Skyscraper case with 64 subdo- Fig. 4.10. Skyscraper case with 128 subdo- 

mains. 16 Ritz vectors spanning the coarse space, mains. 16 Ritz vectors spanning the coarse space. 



In Table |4.4| we report the maximal residual of Ritz pairs for both cases that are 
extracted from the Krylov subspace when solving the system preconditioned only by 
RAS. 

Note that the projection matrix E is large in the case that the decomposition 
has 64 or 128 subdomains. As a consequence, computing LU factorization of E~ x 
is costly and impairs the gains in the number of iterations. Hence, we attempt to 
compute the incomplete LU factorization of E, which is cheaper to compute than the 
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Fig. 4.11. Continuous case with 16 subdo- 



FlG. 4.12. Continuous case with 32 subdo- 



mains. 15 Ritz vectors spanning the coarse space, mains. 15 Ritz vectors spanning the coarse space. 
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Iteration number 
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Iteration number 



Fig. 4.13. Continuous case with 64 subdo- Fig. 4.14. Continuous case with 128 subdo- 

raains. 16 Ritz vectors spanning the coarse space, mains. 15 Ritz vectors spanning the coarse space. 

Table 4.4 
Maximal residual of Ritz pairs for the skyscraper and continuous cases. 



Nparts 


16 


32 


64 


128 


skyscraper 


2.564719e-09 


9.162259e-08 


1.577898c-08 


4.747077e-09 


continuous 


4.439121e-03 


5.242112e-03 


4.749042e-03 


1.326165e-03 



LU factorization. Assume L and U are factors of an incomplete LU factorization with 
no fill-in (ILU(0)) of E. Note that E has a sparse structure because of the sparse 
structure of Z . So L and U are sparse as well. This means that it is very cheap to 
solve {LU)x — y. In this way, E is actually replaced by LU. 

From Figure |4.15| and Figure |4.16[ it appears that the perturbation in E has 
almost no impact on P4 and Pq for the skycraper case, but has severe impact on Pry 
that leads to a stagnation in the convergence. Figure [4. 17| and Figure |4~18 show that 
Pa and Pc are also stable for the continuous case. As is shown in Figure 4.18[ Pn is 
still unstable when the perturbation in E is not small enough (see the second column 



in Table 4.6 1. However, Figure 4.17 shows that Pd is stable, since the perturbation 



is small enough such that LU is almost same as E (see the first column in Table 4.6) 
Table [4~5] and Table [476] present the distance between LU and E for both cases 
respectively. 
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Fig. 4.15. Skyscraper case with 64 subdo- 
mains. E is replaced by LU. 



Fig. 4.16. Skyscraper case with 128 subdo- 
raains. E is replaced by LU . 



Table 4.5 
The distance between LU and E for the skyscraper case. 



Nparts 


64 


128 


\\E(LU)- L - 1\\ 2 


3.8588e-08 


8.9712e+02 


\\(LU)-"E-I\\ 2 


4.1433e-10 


8.6292 
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Fig. 4.17. Continuous case with 64 subdo- 
raains. E is replaced by LU . 



Fig. 4.18. Continuous case with 128 sub- 
domains. E is replaced by LU. 



Table 4.6 
The distance between LU and E for the continuous case. 



Nparts 


64 


128 


\\E{LU)- l -I\\ 2 


4.1008e-14 


5.4090e-01 


\\{LU)- l E-I\\ 2 


1.3890e-15 


1.0215e-01 



5. Conclusion. We presented a perturbation analysis on the deflation, coarse 
correction and adapted deflation preconditioned when the inexact coarse space and 
inverse of projection matrix are applied for the construction of the preconditioned. 
Our analysis shows that in exact arithmetic the spectrum of the system preconditioned 
by the three prcconditioncrs is impacted by the angle between the exact coarse space 
and the perturbed one. Moreover, we prove that the coarse correction and adapted 
deflation preconditioners are insensitive to the perturbation of the projection matrix, 
whereas the deflation preconditioner is sensitive. Numerical results of the different 
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test cases confirm the perturbation analysis. 
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